remove(list=ls())
set.seed(8675309)

library(ggplot2)
library(dplyr)
library(lubridate)
library(stargazer)
library(gridExtra)
library(scales)
library(countrycode)
library(rio)
library(stargazer)
library(tidyr)
library(knitr)
library(psych)
library(kableExtra)
library(xtable)
library(lubridate)
library(scales)
library(ggpubr)



# Import survey data -----------------------------------------------------
# set to your local WD
setwd("C:/Users/Mansfield/Dropbox/Research/China and IOs/Data and Analysis/Replication")

brazil <- read.csv("IO Executives- Reputation (Portuguese)_December 25, 2023_20.33.csv", stringsAsFactors = FALSE)
brazil$survey <- "Brazil"
france <- read.csv("IO Executives- Reputation (French)_December 25, 2023_20.36.csv", stringsAsFactors = FALSE)
france$survey <- "France"

data <- rbind(brazil, france)

# Completion Rates -----------------------------------------------------
#completion rate
completed_100_count <- data %>%
  filter(Progress == 100) %>%
  summarise(count = n())

#attention check
count_red_green <- data %>%
  filter(screener == "Red,Green") %>%  
  group_by(survey)  %>% 
  summarise(count = n())

# Filter Respondents -----------------------------------------------------

#remove those that didn't complete the survey 
data <- data %>%
  filter(Progress == 100)

#remove those that didn't pass the attention check
dt <- data %>%
  filter(screener == "Red,Green") 

dt %>% 
  group_by(survey) %>% 
  summarise(n())

# Dependent Variable: Countries -----------------------------------------------------
#reputation: How does leading change reputation? 
dt$reputation_num[dt$Reputation1_1 == "Very positive effect"] = 5
dt$reputation_num[dt$Reputation1_1 == "Somewhat positive effect"] = 4
dt$reputation_num[dt$Reputation1_1 == "Neither negative nor positive effect"] = 3
dt$reputation_num[dt$Reputation1_1 == "Somewhat negative effect"] = 2
dt$reputation_num[dt$Reputation1_1 == "Very negative effect"] = 1

dt$reputation_un[dt$Reputation1_2 == "Very positive effect"] = 5
dt$reputation_un[dt$Reputation1_2 == "Somewhat positive effect"] = 4
dt$reputation_un[dt$Reputation1_2 == "Neither negative nor positive effect"] = 3
dt$reputation_un[dt$Reputation1_2 == "Somewhat negative effect"] = 2
dt$reputation_un[dt$Reputation1_2 == "Very negative effect"] = 1

#reputation: How much do you approve after UN election? 
dt$reputation2_num[dt$Reputation2_1 == "Strongly approve"] = 5
dt$reputation2_num[dt$Reputation2_1 == "Somewhat approve"] = 4
dt$reputation2_num[dt$Reputation2_1 == "Neither approve nor disapprove"] = 3
dt$reputation2_num[dt$Reputation2_1 == "Somewhat disapprove"] = 2
dt$reputation2_num[dt$Reputation2_1 == "Strongly disapprove"] = 1

dt$reputation_un2[dt$Reputation2_2 == "Strongly approve"] = 5
dt$reputation_un2[dt$Reputation2_2 == "Somewhat approve"] = 4
dt$reputation_un2[dt$Reputation2_2 == "Neither approve nor disapprove"] = 3

dt$reputation_un2[dt$Reputation2_2 == "Somewhat disapprove"] = 2
dt$reputation_un2[dt$Reputation2_2 == "Strongly disapprove"] = 1

#On a scale of 1 (no confidence) to 5 (full confidence) how much confidence do you have in each of:
dt$legitimacy_chn[dt$Legitimacy_1 == "Very confident"] = 5
dt$legitimacy_chn[dt$Legitimacy_1 == "Somewhat confident"] = 4
dt$legitimacy_chn[dt$Legitimacy_1 == "Neither confident nor unconfident"] = 3
dt$legitimacy_chn[dt$Legitimacy_1 == "Not very confident"] = 2
dt$legitimacy_chn[dt$Legitimacy_1 == "No confidence"] = 1

dt$legitimacy_usa[dt$Legitimacy_2 == "Very confident"] = 5
dt$legitimacy_usa[dt$Legitimacy_2 == "Somewhat confident"] = 4
dt$legitimacy_usa[dt$Legitimacy_2 == "Neither confident nor unconfident"] = 3
dt$legitimacy_usa[dt$Legitimacy_2 == "Not very confident"] = 2
dt$legitimacy_usa[dt$Legitimacy_2 == "No confidence"] = 1

dt$legitimacy_su[dt$Legitimacy_3 == "Very confident"] = 5
dt$legitimacy_su[dt$Legitimacy_3 == "Somewhat confident"] = 4
dt$legitimacy_su[dt$Legitimacy_3 == "Neither confident nor unconfident"] = 3
dt$legitimacy_su[dt$Legitimacy_3 == "Not very confident"] = 2
dt$legitimacy_su[dt$Legitimacy_3 == "No confidence"] = 1

dt$legitimacy_un[dt$Legitimacy_4 == "Very confident"] = 5
dt$legitimacy_un[dt$Legitimacy_4 == "Somewhat confident"] = 4
dt$legitimacy_un[dt$Legitimacy_4 == "Neither confident nor unconfident"] = 3
dt$legitimacy_un[dt$Legitimacy_4 == "Not very confident"] = 2
dt$legitimacy_un[dt$Legitimacy_4 == "No confidence"] = 1

#trust: For each of the following, how much do you tend to trust it or tend not to trust it?
dt$trust_chn[dt$Trust_1 == "Trust completely"] = 5
dt$trust_chn[dt$Trust_1 == "Somewhat trust"] = 4
dt$trust_chn[dt$Trust_1 == "Neither Trust nor distrust"] = 3
dt$trust_chn[dt$Trust_1 == "Mostly distrust"] = 2
dt$trust_chn[dt$Trust_1 == "Do not trust at all"] = 1

dt$trust_usa[dt$Trust_2 == "Trust completely"] = 5
dt$trust_usa[dt$Trust_2 == "Somewhat trust"] = 4
dt$trust_usa[dt$Trust_2 == "Neither Trust nor distrust"] = 3
dt$trust_usa[dt$Trust_2 == "Mostly distrust"] = 2
dt$trust_usa[dt$Trust_2 == "Do not trust at all"] = 1

dt$trust_su[dt$Trust_3 == "Trust completely"] = 5
dt$trust_su[dt$Trust_3 == "Somewhat trust"] = 4
dt$trust_su[dt$Trust_3 == "Neither Trust nor distrust"] = 3
dt$trust_su[dt$Trust_3 == "Mostly distrust"] = 2
dt$trust_su[dt$Trust_3 == "Do not trust at all"] = 1

dt$trust_un[dt$Trust_4 == "Trust completely"] = 5
dt$trust_un[dt$Trust_4 == "Somewhat trust"] = 4
dt$trust_un[dt$Trust_4 == "Neither Trust nor distrust"] = 3
dt$trust_un[dt$Trust_4 == "Mostly distrust"] = 2
dt$trust_un[dt$Trust_4 == "Do not trust at all"] = 1

#leadership: Suppose either China or the United States will be the most powerful nation in the world in ten years. 
dt$leader[dt$Leadership == "Strongly prefer China"] = 5
dt$leader[dt$Leadership == "Somewhat prefer China"] = 4
dt$leader[dt$Leadership == "Prefer neither China nor the United States"] = 3
dt$leader[dt$Leadership == "Somewhat prefer the United States"] = 2
dt$leader[dt$Leadership == "Strongly prefer the United States"] = 1

#cooperation: this entity poses an opportunity for cooperation
dt$cooperation_chn[dt$Cooperation_1 == "Strongly agree"] = 5
dt$cooperation_chn[dt$Cooperation_1 == "Somewhat agree"] = 4
dt$cooperation_chn[dt$Cooperation_1 == "Neither agree nor disagree"] = 3
dt$cooperation_chn[dt$Cooperation_1 == "Somewhat disagree"] = 2
dt$cooperation_chn[dt$Cooperation_1 == "Strongly disagree"] = 1

dt$cooperation_usa[dt$Cooperation_2 == "Strongly agree"] = 5
dt$cooperation_usa[dt$Cooperation_2 == "Somewhat agree"] = 4
dt$cooperation_usa[dt$Cooperation_2 == "Neither agree nor disagree"] = 3
dt$cooperation_usa[dt$Cooperation_2 == "Somewhat disagree"] = 2
dt$cooperation_usa[dt$Cooperation_2 == "Strongly disagree"] = 1

dt$cooperation_su[dt$Cooperation_3 == "Strongly agree"] = 5
dt$cooperation_su[dt$Cooperation_3 == "Somewhat agree"] = 4
dt$cooperation_su[dt$Cooperation_3 == "Neither agree nor disagree"] = 3
dt$cooperation_su[dt$Cooperation_3 == "Somewhat disagree"] = 2
dt$cooperation_su[dt$Cooperation_3 == "Strongly disagree"] = 1

#policy instruments: diplomacy 
dt$diplomacy_usa[dt$Instruments_1 == "Totally acceptable"] = 5
dt$diplomacy_usa[dt$Instruments_1 == "Somewhat acceptable"] = 4
dt$diplomacy_usa[dt$Instruments_1 == "Neither acceptable nor unacceptable"] = 3
dt$diplomacy_usa[dt$Instruments_1 == "Somewhat unacceptable"] = 2
dt$diplomacy_usa[dt$Instruments_1 == "Totally unacceptable"] = 1

dt$diplomacy_chn[dt$Instruments_2 == "Totally acceptable"] = 5
dt$diplomacy_chn[dt$Instruments_2 == "Somewhat acceptable"] = 4
dt$diplomacy_chn[dt$Instruments_2 == "Neither acceptable nor unacceptable"] = 3
dt$diplomacy_chn[dt$Instruments_2 == "Somewhat unacceptable"] = 2
dt$diplomacy_chn[dt$Instruments_2 == "Totally unacceptable"] = 1

#policy instruments: receive aid and infrastructure 
dt$aid_chn[dt$Instruments_3 == "Totally acceptable"] = 5
dt$aid_chn[dt$Instruments_3 == "Somewhat acceptable"] = 4
dt$aid_chn[dt$Instruments_3 == "Neither acceptable nor unacceptable"] = 3
dt$aid_chn[dt$Instruments_3 == "Somewhat unacceptable"] = 2
dt$aid_chn[dt$Instruments_3 == "Totally unacceptable"] = 1

dt$aid_usa[dt$Instruments_4 == "Totally acceptable"] = 5
dt$aid_usa[dt$Instruments_4 == "Somewhat acceptable"] = 4
dt$aid_usa[dt$Instruments_4 == "Neither acceptable nor unacceptable"] = 3
dt$aid_usa[dt$Instruments_4 == "Somewhat unacceptable"] = 2
dt$aid_usa[dt$Instruments_4 == "Totally unacceptable"] = 1

#policy instruments: business partnerships
dt$business_chn[dt$Instruments_5 == "Totally acceptable"] = 5
dt$business_chn[dt$Instruments_5 == "Somewhat acceptable"] = 4
dt$business_chn[dt$Instruments_5 == "Neither acceptable nor unacceptable"] = 3
dt$business_chn[dt$Instruments_5 == "Somewhat unacceptable"] = 2
dt$business_chn[dt$Instruments_5 == "Totally unacceptable"] = 1

dt$business_usa[dt$Instruments_6 == "Totally acceptable"] = 5
dt$business_usa[dt$Instruments_6 == "Somewhat acceptable"] = 4
dt$business_usa[dt$Instruments_6 == "Neither acceptable nor unacceptable"] = 3
dt$business_usa[dt$Instruments_6 == "Somewhat unacceptable"] = 2
dt$business_usa[dt$Instruments_6 == "Totally unacceptable"] = 1

# Create Reputation Index --------------------------------------------
#Kronbach's alpha
test1ch <- dt %>% dplyr::select(legitimacy_chn, trust_chn)
test1us <- dt %>% dplyr::select(legitimacy_usa, trust_usa)
test1su <- dt %>% dplyr::select(legitimacy_su, trust_su)
psych::alpha(test1ch)
psych::alpha(test1us)
psych::alpha(test1su)

test2us <- dt %>% dplyr::select(business_usa, aid_usa, diplomacy_usa, cooperation_usa)
test2chn<- dt %>% dplyr::select(business_chn, aid_chn, diplomacy_chn, cooperation_chn)
psych::alpha(test2us)
psych::alpha(test2chn)

test3 <- dt %>% dplyr::select(reputation_num, reputation2_num)
psych::alpha(test3)

test4_un <- dt %>% dplyr::select(reputation_un, reputation_un2, legitimacy_un, trust_un)
psych::alpha(test4_un) 

#calculate average across DV
dt$image_chn <-  rowMeans(dt[ ,c("legitimacy_chn", "trust_chn")], na.rm = TRUE)
dt$fp_cooperation_chn <-  rowMeans(dt[ ,c("business_chn", "aid_chn", "diplomacy_chn", "cooperation_chn")], na.rm = TRUE)

dt$image_usa <-  rowMeans(dt[ ,c("legitimacy_usa", "trust_usa")], na.rm = TRUE)
dt$fp_cooperation_usa <-  rowMeans(dt[ ,c("business_usa", "aid_usa", "diplomacy_usa", "cooperation_usa")], na.rm = TRUE)

dt$image_su <-  rowMeans(dt[ ,c("legitimacy_su", "trust_su")], na.rm = TRUE)

dt$io_image <-  rowMeans(dt[ ,c("reputation_un", "reputation_un2", 
                                "legitimacy_un", "trust_un")], na.rm = TRUE)
#rescale function
rescale <- function(x){
  return((x-min(x,na.rm=TRUE))/(max(x-min(x,na.rm=TRUE),na.rm=TRUE)))
}

# Scale average measures for easy interpretation (0-100)
dt$image_chn_index <- rescale(dt$image_chn)
dt$fp_cooperation_chn_index <- rescale(dt$fp_cooperation_chn)

dt$image_usa_index <- rescale(dt$image_usa)
dt$fp_cooperation_usa_index <- rescale(dt$fp_cooperation_usa)

dt$image_su_index <- rescale(dt$image_su)

dt$un_index <- rescale(dt$io_image)

#scale individual variables
dt$leader <- rescale(dt$leader)
dt$rep1 <- rescale(dt$reputation_num)
dt$rep2 <- rescale(dt$reputation2_num)

#scale individual variables: CHN
dt$leg <- rescale(dt$legitimacy_chn)
dt$trust <- rescale(dt$trust_chn)
dt$coop <- rescale(dt$cooperation_chn)
dt$diplomacy <- rescale(dt$diplomacy_chn)
dt$aid <- rescale(dt$aid_chn)
dt$biz <- rescale(dt$business_chn)

#scale individual variables: USA
dt$leg_usa <- rescale(dt$legitimacy_usa)
dt$trust_usa <- rescale(dt$trust_usa)
dt$coop_usa <- rescale(dt$cooperation_usa)
dt$diplomacy_usa<- rescale(dt$diplomacy_usa)
dt$aid_usa <- rescale(dt$aid_usa)
dt$biz_usa <- rescale(dt$business_usa)

#scale individual variables: Swiss
dt$leg_su <- rescale(dt$legitimacy_su)
dt$trust_su <- rescale(dt$trust_su)
dt$coop_su <- rescale(dt$cooperation_su)

#scale individual variables: UN
dt$leg_un <- rescale(dt$legitimacy_un)
dt$trust_un <- rescale(dt$trust_un)
dt$rep1_un <- rescale(dt$reputation_un)
dt$rep2_un <- rescale(dt$reputation_un2)

# Control Variables -----------------------------------------------------
#education
dt$education[dt$EDUCATION == "Post-graduate degree"] = 5
dt$education[dt$EDUCATION == "College/university graduate"] = 4
dt$education[dt$EDUCATION == "Some college/Associate’s degree"] = 3
dt$education[dt$EDUCATION == "High school graduate/GED"] = 2
dt$education[dt$EDUCATION == "Elementary or some high school"] = 1

#income
dt$income[dt$Income == "$150,000 or more"] = 6
dt$income[dt$Income == "$100,000-$149,999"] = 5
dt$income[dt$Income == "$75,000-$99,000"] = 4
dt$income[dt$Income == "$50,000-$74,999"] = 3
dt$income[dt$Income == "$25,000-$49,999"] = 2
dt$income[dt$Income == "Less than $25,000"] = 1
dt$income[dt$Income == "Prefer not to say"] = NA

#China enemy 
dt$chn_friendly[dt$Frenemy_1 == "Enemy"] = 5
dt$chn_friendly[dt$Frenemy_1 == "Unfriendly"] = 4
dt$chn_friendly[dt$Frenemy_1 == "Not sure"] = 3
dt$chn_friendly[dt$Frenemy_1 == "Friendly"] = 2
dt$chn_friendly[dt$Frenemy_1 == "Ally"] = 1

#China threat
dt$chn_threat[dt$Threat.Perception_1 == "Definitely agree"] = 5
dt$chn_threat[dt$Threat.Perception_1 == "Somewhat agree"] = 4
dt$chn_threat[dt$Threat.Perception_1 == "Neither agree nor disagree"] = 3
dt$chn_threat[dt$Threat.Perception_1 == "Somewhat disagree"] = 2
dt$chn_threat[dt$Threat.Perception_1 == "Definitely disagree"] = 1

#US enemy 
dt$usa_friendly[dt$Frenemy_2 == "Enemy"] = 5
dt$usa_friendly[dt$Frenemy_2 == "Unfriendly"] = 4
dt$usa_friendly[dt$Frenemy_2 == "Not sure"] = 3
dt$usa_friendly[dt$Frenemy_2 == "Friendly"] = 2
dt$usa_friendly[dt$Frenemy_2 == "Ally"] = 1

#US threat
dt$usa_threat[dt$Threat.Perception_2 == "Definitely agree"] = 5
dt$usa_threat[dt$Threat.Perception_2 == "Somewhat agree"] = 4
dt$usa_threat[dt$Threat.Perception_2 == "Neither agree nor disagree"] = 3
dt$usa_threat[dt$Threat.Perception_2 == "Somewhat disagree"] = 2
dt$usa_threat[dt$Threat.Perception_2 == "Definitely disagree"] = 1

#Trust in government
dt$trustingov[dt$TRUST_GOVT == "Just about always"] = 3
dt$trustingov[dt$TRUST_GOVT == "Most of the time"] = 2
dt$trustingov[dt$TRUST_GOVT == "Only some of the time"] = 1

#Gender
dt$male[dt$GENDER == "Male"] = 2
dt$male[dt$GENDER == "Female"] = 1
dt$male[dt$GENDER == "Neither/Prefer not to say"] = NA

# Ideology
dt$ideol[dt$IDEOLOGY == "Extremely conservative"] = 5
dt$ideol[dt$IDEOLOGY == "Conservative"] = 4
dt$ideol[dt$IDEOLOGY == "Moderate, middle of the road"] = 3
dt$ideol[dt$IDEOLOGY == "Liberal"] = 2
dt$ideol[dt$IDEOLOGY == "Extremely liberal"] = 1

# Political Interest
dt$polinterest[dt$POL_INTEREST == "Most of the time"] = 4
dt$polinterest[dt$POL_INTEREST == "Some of the time"] = 3
dt$polinterest[dt$POL_INTEREST == "Only now and then"] = 2
dt$polinterest[dt$POL_INTEREST == "Hardly at all"] = 1

# Foreign policy orientation
dt$fp1 <- car::recode(dt$FP_ORIENTATION_1,
      "'Definitely agree'=5; 'Somewhat agree'=4; 'Neither agree nor disagree'=3; 'Somewhat disagree'=2; 'Definitely disagree'=1")
dt$fp2 <- car::recode(dt$FP_ORIENTATION_2,
    "'Definitely agree'=5; 'Somewhat agree'=4; 'Neither agree nor disagree'=3; 'Somewhat disagree'=2; 'Definitely disagree'=1")
#reverse coded 
dt$fp3 <- car::recode(dt$FP_ORIENTATION_3,
    "'Definitely agree'=1; 'Somewhat agree'=2; 'Neither agree nor disagree'=3; 'Somewhat disagree'=4; 'Definitely disagree'=5")
#reverse coded
dt$fp4 <-  car::recode(dt$FP_ORIENTATION_4,
    "'Definitely agree'=1; 'Somewhat agree'=2; 'Neither agree nor disagree'=3; 'Somewhat disagree'=4; 'Definitely disagree'=5")
dt$fp_index <-  rowMeans(dt[ ,c("fp1", "fp2", "fp3", "fp4")], na.rm = TRUE)

# Age : remove 3 obs
dt <- dt %>% 
  filter(Age < 83)

# DESCRIPTIVES: CHN Threat -----------------------------------------------------
# Reshaping the data to have both chn_friendly and chn_threat in long format
chn_long <- dt %>%
  pivot_longer(cols = c(chn_friendly, chn_threat), 
               names_to = "variable", values_to = "value")

# Calculating mean, lower.ci, and upper.ci for each survey and variable
chn_summary <- chn_long %>%
  group_by(survey, variable) %>%
  summarize(
    mean = mean(value, na.rm = TRUE),
    lower.ci = mean - qt(0.975, df=n()-1) * sd(value, na.rm = TRUE) / sqrt(n()),
    upper.ci = mean + qt(0.975, df=n()-1) * sd(value, na.rm = TRUE) / sqrt(n())
  )

# Define new professional colors
new_colors <- c("chn_threat" = "#1F77B4",  # Blue for China Threat
                "chn_enemy" = "#7F7F7F")   # Gray for China Enemy

chn_threat_plot <- ggplot(chn_summary, aes(x = survey, y = mean, fill = variable)) + 
  geom_bar(stat = "identity", position = position_dodge(width = 0.9)) +
  labs(x = "Country", y = "China Threat Perception (1-5)", 
       fill = "Variable",
       title = "Brazil and France: China Threat Perception") +
  scale_fill_manual(values = new_colors,
                    labels = c("CHN Threat", "CHN Enemy")) +  # Adjusted order to match variable names
  theme_bw() +
  coord_cartesian(ylim = c(1, 3.5)) +
  theme(
    plot.caption = element_text(hjust = 0, family = "Times New Roman"),
    legend.position = "bottom",
    legend.title = element_text(family = "Times New Roman"),
    legend.text = element_text(family = "Times New Roman"),
    axis.title = element_text(family = "Times New Roman"),
    axis.text = element_text(family = "Times New Roman", size = 12)
  )

# Figure 1
chn_threat_plot

# TREATMENTS -----------------------------------------------------
dt$TREATMENT2[dt$TREATMENT == 1] = "CHN"
dt$TREATMENT2[dt$TREATMENT == 2] = "USA"
dt$TREATMENT2[dt$TREATMENT == "Control"] = "Swiss"

# DESCRIPTIVES: Balance -----------------------------------------------------
table(dt$TREATMENT)
table(dt$TREATMENT2)

dt_chn_france <- dt %>% filter(survey == "France")
dt_chn_brazil <- dt %>% filter(survey == "Brazil")

table(dt_chn_france$TREATMENT2) #Distribution of treatment to check RA
table(dt_chn_france$TREATMENT2)

vars <- c('Age', 'male',
          'education', 'trustingov', 'polinterest',
          'fp_index', 'ideol', 'income',
          'chn_friendly', "chn_threat", 
          'usa_friendly', "usa_threat")

var_labels <- c("Age", "Male", 
                "Education", "Trust in Government", "Political Interest",
                "Foreign Policy Orientation", "Conservatism", "Income",
                "China Enemy", "China Threat",
                "USA Enemy", "USA Threat")

treats <- c("China", "USA")
treat_conditions <- c(1:2)

balmat <- data.frame(matrix(NA,length(vars)*length(treats), 5))
colnames(balmat) <- c("Var.",  "Treatment", "T-Test P val.", 
                      "Ctrl. Mean", "Treatment Mean") 

j <- c()
for(i in 1:length(vars)){
  a <- rep(var_labels[i], 2)
  j <- append(x = j, values = a)}

balmat[, 1] <- j
balmat[, 2] <- rep(treats, 12)

counter <- 0
for (i in 1:length(vars)){
  for (k in 1:length(treat_conditions)){
    string <- paste('t_test <- t.test(dt_chn_france$', 
                    vars[i],
                    '[dt_chn_france$TREATMENT==', 
                    treat_conditions[k],
                    '], dt_chn_france$', 
                    vars[i],
                    '[dt_chn_france$TREATMENT=="Control"])', 
                    sep = "", collapse = "")
    eval(parse(text=string))
    balmat[(counter+k),3] <- round(t_test$p.value, digits = 3)
    balmat[(counter+k),4] <- round(t_test$estimate[2], digits = 3)
    balmat[(counter+k),5] <- round(t_test$estimate[1], digits = 3)
    
  }
  counter <- counter+2
}

balmat[balmat$t_pval<.1, ] 

# Table A3
xtable(balmat, 
       font.size = "tiny", caption = "Balance Tests, France Sample")


balmat <- data.frame(matrix(NA,length(vars)*length(treats), 5))
colnames(balmat) <- c("Var.",  "Treatment", "T-Test P val.", 
                      "Ctrl. Mean", "Treatment Mean") 

j <- c()
for(i in 1:length(vars)){
  a <- rep(var_labels[i], 2)
  j <- append(x = j, values = a)}

balmat[, 1] <- j
balmat[, 2] <- rep(treats, 12)

counter <- 0
for (i in 1:length(vars)){
  for (k in 1:length(treat_conditions)){
    string <- paste('t_test <- t.test(dt_chn_brazil$', 
                    vars[i],
                    '[dt_chn_brazil$TREATMENT==', 
                    treat_conditions[k],
                    '], dt_chn_brazil$', 
                    vars[i],
                    '[dt_chn_brazil$TREATMENT=="Control"])', 
                    sep = "", collapse = "")
    eval(parse(text=string))
    balmat[(counter+k),3] <- round(t_test$p.value, digits = 3)
    balmat[(counter+k),4] <- round(t_test$estimate[2], digits = 3)
    balmat[(counter+k),5] <- round(t_test$estimate[1], digits = 3)
    
  }
  counter <- counter+2
}

balmat[balmat$t_pval<.1, ] 

# Table A-4
xtable(balmat, 
       font.size = "tiny", caption = "Balance Tests, Brazil Sample")

# DESCRIPTIVES: Summary stats -----------------------------------------------------
sum_stats <- data.frame(matrix(NA,length(vars), 8))
colnames(sum_stats) <- c("Var.", "Min.", "1st Qu.",  "Median",
                         "Mean", "3rd Qu.",    "Max.", "Num. Missing" )
sum_stats[, 1] <- var_labels

for (i in 1:length(vars)){
  string <- paste('sum <- summary(dt_chn_france$',vars[i], ')',
                  sep = "", collapse = "")
  eval(parse(text=string))
  sum_stats[i, 2:7] <- sum
  
}

missing_vars_fra <- dt_chn_france %>% 
  dplyr::select(Age, male, education,
                trustingov, polinterest,
                fp_index, ideol, income,
                chn_friendly, chn_threat,
                usa_friendly, usa_threat)


missing_fra <- sapply(missing_vars_fra, FUN = function(x) sum(is.na(x)))
sum_stats[, 8] <- missing_fra

# Table A-1
xtable(sum_stats, caption = "Summary Statistics, France Sample", 
       digits = c(0, 0, 0, 2, 2, 2, 2, 0, 0))

sum_stats <- data.frame(matrix(NA,length(vars), 8))
colnames(sum_stats) <- c("Var.", "Min.", "1st Qu.",  "Median",
                         "Mean", "3rd Qu.",    "Max.", "Num. Missing" )
sum_stats[, 1] <- var_labels

for (i in 1:length(vars)){
  string <- paste('sum <- summary(dt_chn_brazil$',vars[i], ')',
                  sep = "", collapse = "")
  eval(parse(text=string))
  sum_stats[i, 2:7] <- sum
  
}

missing_vars_bra <- dt_chn_brazil %>% 
  dplyr::select(Age, male, education,
                trustingov, polinterest,
                fp_index, ideol, income,
                chn_friendly, chn_threat,
                usa_friendly, usa_threat)

missing_bra <- sapply(missing_vars_bra, FUN = function(x) sum(is.na(x)))
sum_stats[, 8] <- missing_bra

# Table A-2
xtable(sum_stats, caption = "Summary Statistics, Brazil Sample", 
       digits = c(0, 0, 0, 2, 2, 2, 2, 0, 0))

# Model: CHN Overall -----------------------------------------------------
dt_chn_france$TREATMENT_relevel2 <- relevel(as.factor(dt_chn_france$TREATMENT), 
                                            ref = "Control")
dt_chn_brazil$TREATMENT_relevel2 <- relevel(as.factor(dt_chn_brazil$TREATMENT), 
                                            ref = "Control")
#France and Brazil: Index
mod1 <- lm(image_chn_index ~ 
             TREATMENT_relevel2 + 
             education + 
             income + 
             Age + male +
             trustingov + polinterest +
             fp_index + ideol +
             chn_friendly + chn_threat +
             usa_friendly + usa_threat, 
           data = dt_chn_france)
mod1c_f <- lm(image_chn_index ~ 
             TREATMENT_relevel2 + 
             education + 
             income + 
             Age + male +
             trustingov + polinterest +
             fp_index + ideol +
             chn_friendly + chn_threat +
             usa_friendly + usa_threat, 
           data = dt_chn_france)

mod2 <- lm(fp_cooperation_chn_index ~ 
             TREATMENT_relevel2 + 
             education + 
             income + 
             Age + male +
             trustingov + polinterest +
             fp_index + ideol +
             chn_friendly + chn_threat +
             usa_friendly + usa_threat, 
           data = dt_chn_france)

# Table A-5
stargazer(mod1, mod2, title = "China Image (France)",
          column.labels = c("Image Index", "FP Cooperation Index"),
          covariate.labels = c("China", "US", "Education", 
                               "Income", "Age", "Male", "Trust in Govt", "Pol Interest", 
                               "FP Orientation", "Conservatism",  
                               "China Enemy", "China Threat", "USA Enemy", "USA Threat"),
          star.cutoffs = c(.05, .01, .001))

mod3 <- lm(image_chn_index ~ 
             TREATMENT_relevel2 + 
             education + 
             income + 
             Age + male +
             trustingov + polinterest +
             fp_index + ideol +
             chn_friendly + chn_threat +
             usa_friendly + usa_threat, 
           data = dt_chn_brazil)
mod3c_b <- lm(image_chn_index ~ 
             TREATMENT_relevel2 + 
             education + 
             income + 
             Age + male +
             trustingov + polinterest +
             fp_index + ideol +
             chn_friendly + chn_threat +
             usa_friendly + usa_threat, 
           data = dt_chn_brazil)

mod4 <- lm(fp_cooperation_chn_index ~ 
             TREATMENT_relevel2 + 
             education + 
             income + 
             Age + male +
             trustingov + polinterest +
             fp_index + ideol +
             chn_friendly + chn_threat +
             usa_friendly + usa_threat, 
           data = dt_chn_brazil)

# Table A-6
stargazer(mod3, mod4, title = "China Image (Brazil)",
          column.labels = c("Image Index", "FP Cooperation Index"),
          covariate.labels = c("China", "US", "Education", 
                               "Income", "Age", "Male", "Trust in Govt", "Pol Interest", 
                               "FP Orientation", "Conservatism",  
                               "China Enemy", "China Threat", "USA Enemy", "USA Threat"),
          star.cutoffs = c(.05, .01, .001))

p_values1 <- c(summary(mod3)$coefficients["TREATMENT_relevel21", "Pr(>|t|)"])

#France: All DVs
m1 <- lm(leg ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)
m2 <- lm(trust ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)
m3 <- lm(coop ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)
m4 <- lm(diplomacy ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)
m5 <- lm(aid ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)
m6 <- lm(biz ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)
m7 <- lm(leader ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)

# Table A-11
stargazer(m1, m2, m3, m4, m5, m6, m7, title = "Experimental Results: France (China UN)",
          column.labels = c("Legitimacy", "Trust", "Cooperation",
                            "Diplomacy", "Aid", "Business", "Leader"),
          covariate.labels = c("China", "US", "Education", 
          "Income", "Age", "Male", "Trust in Govt", "Pol Interest", 
          "FP Orientation", "Conservatism",  
          "China Enemy", "China Threat", "USA Enemy", "USA Threat"),
          star.cutoffs = c(.05, .01, .001))

#Brazil: All DVs
m1 <- lm(leg ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)
m2 <- lm(trust ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)
m3 <- lm(coop ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)
m4 <- lm(diplomacy ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)
m5 <- lm(aid ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)
m6 <- lm(biz ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)
m7 <- lm(leader ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)

# Table A-12
stargazer(m1, m2, m3, m4, m5, m6, m7, title = "Experimental Results: Brazil (China UN)",
          column.labels = c("Legitimacy", "Trust", "Cooperation", "Diplomacy", "Aid", "Business", "Leader"),
          covariate.labels = c("China", "US", "Education", 
                               "Income", "Age", "Male", "Trust in Govt", "Pol Interest", 
                               "FP Orientation", "Conservatism",  
                               "China Enemy", "China Threat", "USA Enemy", "USA Threat"),
          star.cutoffs = c(.05, .01, .001))

# Model: USA Overall -----------------------------------------------------
#France and Brazil: Index
mod1 <- lm(image_usa_index ~ TREATMENT_relevel2 + 
             education + 
             income + 
             Age + male +
             trustingov + polinterest +
             fp_index + ideol +
             chn_friendly + chn_threat +
             usa_friendly + usa_threat, data = dt_chn_france)
mod1f <- lm(image_usa_index ~ TREATMENT_relevel2 + 
             education + 
             income + 
             Age + male +
             trustingov + polinterest +
             fp_index + ideol +
             chn_friendly + chn_threat +
             usa_friendly + usa_threat, data = dt_chn_france)
mod2 <- lm(fp_cooperation_usa_index ~ TREATMENT_relevel2 + 
             education + 
             income + 
             Age + male +
             trustingov + polinterest +
             fp_index + ideol +
             chn_friendly + chn_threat +
             usa_friendly + usa_threat, data = dt_chn_france)

# Table A-7
stargazer(mod1, mod2, title = "USA Image (France)",
          column.labels = c("Image Index", "FP Cooperation Index"),
          covariate.labels = c("China", "US", "Education", 
                               "Income", "Age", "Male", "Trust in Govt", "Pol Interest", 
                               "FP Orientation", "Conservatism",  
                               "China Enemy", "China Threat", "USA Enemy", "USA Threat"),
          star.cutoffs = c(.05, .01, .001))

mod3 <- lm(image_usa_index ~ TREATMENT_relevel2 + 
             education + 
             income + 
             Age + male +
             trustingov + polinterest +
             fp_index + ideol +
             chn_friendly + chn_threat +
             usa_friendly + usa_threat, data = dt_chn_brazil)
mod3b <- lm(image_usa_index ~ TREATMENT_relevel2 + 
             education + 
             income + 
             Age + male +
             trustingov + polinterest +
             fp_index + ideol +
             chn_friendly + chn_threat +
             usa_friendly + usa_threat, data = dt_chn_brazil)
mod4 <- lm(fp_cooperation_usa_index ~ TREATMENT_relevel2 + 
             education + 
             income + 
             Age + male +
             trustingov + polinterest +
             fp_index + ideol +
             chn_friendly + chn_threat +
             usa_friendly + usa_threat, data = dt_chn_brazil)

# Table A-8
stargazer(mod3, mod4, title = "USA Image (Brazil)",
          column.labels = c("Image Index", "FP Cooperation Index"),
          covariate.labels = c("China", "US", "Education", 
                               "Income", "Age", "Male", "Trust in Govt", "Pol Interest", 
                               "FP Orientation", "Conservatism",  
                               "China Enemy", "China Threat", "USA Enemy", "USA Threat"),
          star.cutoffs = c(.05, .01, .001))

#France: All DVs
m1 <- lm(leg_usa ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)
m2 <- lm(trust_usa ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)
m3 <- lm(coop_usa ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)
m4 <- lm(diplomacy_usa ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)
m5 <- lm(aid_usa ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)
m6 <- lm(biz_usa ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)

# Table A-13
stargazer(m1, m2, m3, m4, m5, m6, title = "USA Image (France)",
          column.labels = c("Legitimacy", "Trust", "Cooperation", "Diplomacy", "Aid", "Business", "Leader"),
          covariate.labels = c("China", "US", "Education", 
                               "Income", "Age", "Male", "Trust in Govt", "Pol Interest", 
                               "FP Orientation", "Conservatism",  
                               "China Enemy", "China Threat", "USA Enemy", "USA Threat"),
          star.cutoffs = c(.05, .01, .001))

#Brazil: All DVs
m1 <- lm(leg_usa ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)
m2 <- lm(trust_usa ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)
m3 <- lm(coop_usa ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)
m4 <- lm(diplomacy_usa ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)
m5 <- lm(aid_usa ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)
m6 <- lm(biz_usa ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)

# Table A-14
stargazer(m1, m2, m3, m4, m5, m6, title = "USA Image (Brazil)",
          column.labels = c("Legitimacy", "Trust", "Cooperation", "Diplomacy", "Aid", "Business", "Leader"),
          covariate.labels = c("China", "US", "Education", 
                               "Income", "Age", "Male", "Trust in Govt", "Pol Interest", 
                               "FP Orientation", "Conservatism",  
                               "China Enemy", "China Threat", "USA Enemy", "USA Threat"),
          star.cutoffs = c(.05, .01, .001))

# in-text comparisons
library(lmtest)
library(sandwich)

# compare china image to usa image in france 
coefs <- coef(mod1c_f)["TREATMENT_relevel21"] - coef(mod1f)["TREATMENT_relevel21"]
se <- sqrt(vcovHC(mod1c_f, type = "HC1")["TREATMENT_relevel21", "TREATMENT_relevel21"] +
             vcovHC(mod1f, type = "HC1")["TREATMENT_relevel21", "TREATMENT_relevel21"])

# Compute test statistic
z_score <- coefs / se
p_value <- 2 * (1 - pnorm(abs(z_score)))

# compare china image to usa image in brazil 
coefs <- coef(mod3c_b)["TREATMENT_relevel21"] - coef(mod3b)["TREATMENT_relevel21"]
se <- sqrt(vcovHC(mod3c_b, type = "HC1")["TREATMENT_relevel21", "TREATMENT_relevel21"] +
             vcovHC(mod3b, type = "HC1")["TREATMENT_relevel21", "TREATMENT_relevel21"])

z_score <- coefs / se
p_value <- 2 * (1 - pnorm(abs(z_score)))

# compare china image in france to china image in brazil
coefs <- coef(mod1c_f)["TREATMENT_relevel21"] - coef(mod3c_b)["TREATMENT_relevel21"]
se <- sqrt(vcovHC(mod1c_f, type = "HC1")["TREATMENT_relevel21", "TREATMENT_relevel21"] +
             vcovHC(mod3c_b, type = "HC1")["TREATMENT_relevel21", "TREATMENT_relevel21"])

z_score <- coefs / se
p_value <- 2 * (1 - pnorm(abs(z_score)))


# Model: Switzerland Overall------------------------------------------------
#France and Brazil: Index
mod1 <- lm(image_su_index ~ 
             TREATMENT_relevel2 + 
             education + 
             income + 
             Age + male +
             trustingov + polinterest +
             fp_index + ideol +
             chn_friendly + chn_threat +
             usa_friendly + usa_threat, 
           data = dt_chn_france)

mod1 <- lm(image_su_index ~ 
             TREATMENT_relevel2 + 
             education + 
             income + 
             Age + male +
             trustingov + polinterest +
             fp_index + ideol +
             chn_friendly + chn_threat +
             usa_friendly + usa_threat, 
           data = dt_chn_brazil)

# Table A-9
stargazer(mod1, mod2, title = "Switzerland Image",
          column.labels = c("Image Index (Brazil)", "Image Index (France)"),
          covariate.labels = c("China", "US", "Education", 
                               "Income", "Age", "Male", "Trust in Govt", "Pol Interest", 
                               "FP Orientation", "Conservatism",  
                               "China Enemy", "China Threat", "USA Enemy", "USA Threat"),
          star.cutoffs = c(.05, .01, .001))


#France: All DVs
m1 <- lm(leg_su ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)
m2 <- lm(trust_su ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)
m3 <- lm(coop_su ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)

# Table A-15
stargazer(m1, m2, m3, title = "Switzerland Image (France)",
          column.labels = c("Legitimacy", "Trust", "Cooperation"),
          covariate.labels = c("China", "US", "Education", 
                               "Income", "Age", "Male", "Trust in Govt", "Pol Interest", 
                               "FP Orientation", "Conservatism",  
                               "China Enemy", "China Threat", "USA Enemy", "USA Threat"),
          star.cutoffs = c(.05, .01, .001))

#Brazil: All DVs
m1 <- lm(leg_su ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)
m2 <- lm(trust_su ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)
m3 <- lm(coop_su ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)

# Table A-16
stargazer(m1, m2, m3, title = "Switzerland Image (Brazil)",
          column.labels = c("Legitimacy", "Trust", "Cooperation"),
          covariate.labels = c("China", "US", "Education", 
                               "Income", "Age", "Male", "Trust in Govt", "Pol Interest", 
                               "FP Orientation", "Conservatism",  
                               "China Enemy", "China Threat", "USA Enemy", "USA Threat"),
          star.cutoffs = c(.05, .01, .001))


# Model: IO Legitimacy -----------------------------------------------------
dt$TREATMENT_relevel2 <- relevel(as.factor(dt$TREATMENT), ref = "Control")

#Overall: UN Index
un1 <- lm(un_index ~ TREATMENT_relevel2 + 
            education + 
            income + 
            Age + male +
            trustingov + polinterest +
            fp_index + ideol +
            chn_friendly + chn_threat +
            usa_friendly + usa_threat, data = dt)
un2 <- lm(un_index ~ TREATMENT_relevel2 + 
            education + 
            income + 
            Age + male +
            trustingov + polinterest +
            fp_index + ideol +
            chn_friendly + chn_threat +
            usa_friendly + usa_threat, data = dt_chn_brazil)
un3 <- lm(un_index ~ TREATMENT_relevel2 + 
            education + 
            income + 
            Age + male +
            trustingov + polinterest +
            fp_index + ideol +
            chn_friendly + chn_threat +
            usa_friendly + usa_threat, data = dt_chn_france)

# Table A-10
stargazer(un2, un3, title = "UN Image",
          column.labels = c("UN Index (Brazil)","UN Index (France)"),
          covariate.labels = c("China", "US", "Education", 
                               "Income", "Age", "Male", "Trust in Govt", "Pol Interest", 
                               "FP Orientation", "Conservatism",  
                               "China Enemy", "China Threat", "USA Enemy", "USA Threat"),
          star.cutoffs = c(.05, .01, .001))

# compare UN image after China leadership to US leadership in Brazil
coefs <- coef(un2)["TREATMENT_relevel21"] - coef(un2)["TREATMENT_relevel22"]
se <- sqrt(vcovHC(un2, type = "HC1")["TREATMENT_relevel21", "TREATMENT_relevel21"] +
             vcovHC(un2, type = "HC1")["TREATMENT_relevel22", "TREATMENT_relevel22"] - 
  2 * vcov(un2)["TREATMENT_relevel21", "TREATMENT_relevel22"])

z_score <- coefs / se
p_value <- 2 * (1 - pnorm(abs(z_score)))

# compare UN image after China leadership to US leadership in France
coefs <- coef(un3)["TREATMENT_relevel21"] - coef(un3)["TREATMENT_relevel22"]
se <- sqrt(vcovHC(un3, type = "HC1")["TREATMENT_relevel21", "TREATMENT_relevel21"] +
             vcovHC(un3, type = "HC1")["TREATMENT_relevel22", "TREATMENT_relevel22"] - 
             2 * vcov(un3)["TREATMENT_relevel21", "TREATMENT_relevel22"])

z_score <- coefs / se
p_value <- 2 * (1 - pnorm(abs(z_score)))

# compare UN image in france to UN image in brazil when US leads
coefs <- coef(un3)["TREATMENT_relevel22"] - coef(un2)["TREATMENT_relevel22"]
se <- sqrt(vcovHC(un3, type = "HC1")["TREATMENT_relevel22", "TREATMENT_relevel22"] +
             vcovHC(un2, type = "HC1")["TREATMENT_relevel22", "TREATMENT_relevel22"])

z_score <- coefs / se
p_value <- 2 * (1 - pnorm(abs(z_score)))

#Individual variables -- France
m1 <- lm(leg_un ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)
m2 <- lm(trust_un ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)
m3 <- lm(rep1_un ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)
m4 <- lm(rep2_un ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_france)

# Table A-17
stargazer(m1, m2, m3,m4,  title = "UN Image (France)",
          column.labels = c("Legitimacy", "Trust", "Reputation1","Reputation2"),
          covariate.labels = c("China", "US", "Education", 
                               "Income", "Age", "Male", "Trust in Govt", "Pol Interest", 
                               "FP Orientation", "Conservatism",  
                               "China Enemy", "China Threat", "USA Enemy", "USA Threat"),
          star.cutoffs = c(.05, .01, .001))

# Indivdual variables -- Brazil
m1 <- lm(leg_un ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)
m2 <- lm(trust_un ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)
m3 <- lm(rep1_un ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)
m4 <- lm(rep2_un ~ TREATMENT_relevel2 + 
           education + 
           income + 
           Age + male +
           trustingov + polinterest +
           fp_index + ideol +
           chn_friendly + chn_threat +
           usa_friendly + usa_threat, data = dt_chn_brazil)

# Table A-18
stargazer(m1, m2, m3,m4,  title = "UN Image (Brazil)",
          column.labels = c("Legitimacy", "Trust", "Reputation1","Reputation2"),
          covariate.labels = c("China", "US", "Education", 
                               "Income", "Age", "Male", "Trust in Govt", "Pol Interest", 
                               "FP Orientation", "Conservatism",  
                               "China Enemy", "China Threat", "USA Enemy", "USA Threat"),
          star.cutoffs = c(.05, .01, .001))

# Individual outcomes, country image -----------------------------------------------------
library(modelsummary)

# Define your models
models_france <- list(
  mod1 = lm(leg ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_france),
  mod2 = lm(trust ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_france),
  mod3 = lm(coop ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_france),
  mod4 = lm(diplomacy ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_france),
  mod5 = lm(aid ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_france),
  mod6 = lm(biz ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_france),
  mod7 = lm(leader ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_france)
)

models_france_no_controls <- list(
  mod1 = lm(leg ~ TREATMENT_relevel2, data = dt_chn_france),
  mod2 = lm(trust ~ TREATMENT_relevel2 , data = dt_chn_france),
  mod3 = lm(coop ~ TREATMENT_relevel2, data = dt_chn_france),
  mod4 = lm(diplomacy ~ TREATMENT_relevel2, data = dt_chn_france),
  mod5 = lm(aid ~ TREATMENT_relevel2, data = dt_chn_france),
  mod6 = lm(biz ~ TREATMENT_relevel2, data = dt_chn_france),
  mod7 = lm(leader ~ TREATMENT_relevel2, data = dt_chn_france)
)

models_brazil <- list(
  mod1 = lm(leg ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_brazil),
  mod2 = lm(trust ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_brazil),
  mod3 = lm(coop ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_brazil),
  mod4 = lm(diplomacy ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_brazil),
  mod5 = lm(aid ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_brazil),
  mod6 = lm(biz ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_brazil),
  mod7 = lm(leader ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_brazil)
)

models_brazil_no_controls <- list(
  mod1 = lm(leg ~ TREATMENT_relevel2, data = dt_chn_brazil),
  mod2 = lm(trust ~ TREATMENT_relevel2, data = dt_chn_brazil),
  mod3 = lm(coop ~ TREATMENT_relevel2, data = dt_chn_brazil),
  mod4 = lm(diplomacy ~ TREATMENT_relevel2, data = dt_chn_brazil),
  mod5 = lm(aid ~ TREATMENT_relevel2, data = dt_chn_brazil),
  mod6 = lm(biz ~ TREATMENT_relevel2, data = dt_chn_brazil),
  mod7 = lm(leader ~ TREATMENT_relevel2, data = dt_chn_brazil)
)
# Create a mapping for coefficients
cm <- c('TREATMENT_relevel21' = 'China UN Leadership',
        'TREATMENT_relevel22' = 'US UN Leadership')

cm_no_controls <- c('TREATMENT_relevel21' = 'China UN Leadership',
        'TREATMENT_relevel22' = 'US UN Leadership',
        'education' = 'Education', 
        'income' = 'Income', 
        'Age' = 'Age',
        'male' = 'Male',
        'trustingov' = "Trust Gov",
        'polinterest' = 'Pol Interest',
        'fp_index' = "FP Orientation",
        'ideol' = "Conservatism",
        'chn_friendly' = 'China Enemy',
        'chn_threat' = 'China Threat',
        'usa_friendly' = 'USA Enemy',
        'usa_threat' = 'USA Threat')

# Generate model plot data for France
nb_france <- modelplot(models_france, coef_map = cm, draw = FALSE)
nb_france_no_controls <- modelplot(models_france_no_controls,
                           coef_map = cm_no_controls, draw = F)

# Generate model plot data for Brazil
nb_brazil <- modelplot(models_brazil, coef_map = cm, draw = FALSE)
nb_brazil_no_controls <- modelplot(models_brazil_no_controls,
                                   coef_map = cm_no_controls, draw = F)

# Combine the French and Brazilian data, adding a 'Country' column
nb_combined <- bind_rows(
  mutate(nb_france, Country = "France"),
  mutate(nb_brazil, Country = "Brazil")
)

nb_combined_no_controls <- bind_rows(
  mutate(nb_france_no_controls, Country = "France"),
  mutate(nb_brazil_no_controls, Country = "Brazil")
)

# Assign Colors
new_colors <- c("Confidence" = "#8da0cb", "Trust" = "#fc8d62", 
                          "Cooperation" = "#66c2a5", "Diplomacy" = "#e78ac3", 
                          "Aid" = "#a6d854", "Business" = "#ffd92f", 
                          "Leader" = "#e5c494")


# Renaming the models 
nb_combined <- nb_combined %>% 
  dplyr::mutate(model = dplyr::recode(model, 
                        "mod1" = "Confidence", 
                        "mod2" = "Trust",
                        "mod3" = "Cooperation", 
                        "mod4" = "Diplomacy",
                        "mod5" = "Aid",
                        "mod6" = "Business",
                        "mod7" = "Leader"
  ))

nb_combined_no_controls <- nb_combined_no_controls %>% 
  dplyr::mutate(model = dplyr::recode(model, 
                                      "mod1" = "Confidence", 
                                      "mod2" = "Trust",
                                      "mod3" = "Cooperation", 
                                      "mod4" = "Diplomacy",
                                      "mod5" = "Aid",
                                      "mod6" = "Business",
                                      "mod7" = "Leader"
  ))

nb_combined <- nb_combined %>% 
  dplyr::filter(term == "China UN Leadership" | term == "US UN Leadership")

nb_combined_no_controls <- nb_combined_no_controls %>% 
  dplyr::filter(term == "China UN Leadership" | term == "US UN Leadership")

# Reorder the Country factor to ensure France is first
nb_combined$Country <- factor(nb_combined$Country, levels = c("France", "Brazil"))
nb_combined_no_controls$Country <- factor(nb_combined_no_controls$Country, 
                                          levels = c("France", "Brazil"))


# Create the plot with facet wrap for each country
change_in_coefficients_plot_nonbinary <- ggplot(nb_combined, aes(y = term, x = estimate, 
                                                                 xmin = conf.low, xmax = conf.high, 
                                                                 color = model)) +
  geom_pointrange(position = position_dodge2(width = 0.5, reverse = TRUE), size = 0.25) +
  geom_vline(xintercept = 0) +
  scale_color_manual(values = new_colors) +
  labs(x = "", y = "", 
       title = "",
       caption = "Bars represent 95% confidence intervals for coefficient estimates",
       color = "Model") +
  facet_wrap(~Country, ncol = 1) +
  theme_bw() +
  theme(#plot.title = element_text(hjust = 0.5, family = "Times New Roman"),
        plot.caption = element_text(hjust = 0, family = "Times New Roman"),
        legend.position = "bottom",
        legend.title = element_text(family = "Times New Roman"),
        legend.text = element_text(family = "Times New Roman"),
        axis.title = element_text(family = "Times New Roman"),
        axis.text = element_text(family = "Times New Roman", size = 10))

# Figure A-3
change_in_coefficients_plot_nonbinary

change_in_coefficients_plot_nonbinary_no_controls <- ggplot(nb_combined_no_controls, 
                            aes(y = term, x = estimate, 
                            xmin = conf.low, xmax = conf.high, 
                            color = model)) +
  geom_pointrange(position = position_dodge2(width = 0.5, reverse = TRUE), size = 0.25) +
  geom_vline(xintercept = 0) +
  scale_color_manual(values = new_colors) +
  labs(x = "", y = "", 
       title = "",
       caption = "Bars represent 95% confidence intervals for coefficient estimates",
       color = "Model") +
  facet_wrap(~Country, ncol = 1) +
  theme_bw() +
  theme(#plot.title = element_text(hjust = 0.5, family = "Times New Roman"),
    plot.caption = element_text(hjust = 0, family = "Times New Roman"),
    legend.position = "bottom",
    legend.title = element_text(family = "Times New Roman"),
    legend.text = element_text(family = "Times New Roman"),
    axis.title = element_text(family = "Times New Roman"),
    axis.text = element_text(family = "Times New Roman", size = 10))

# Figure A-6
change_in_coefficients_plot_nonbinary_no_controls

# Save the plot using ggsave for high resolution
# width <- 8 # in inches
# height <- 6 # in inches
# dpi <- 300 # adjust for desired resolution
# ggsave("hypothesis1_updated_RR.png", 
#        plot = change_in_coefficients_plot_nonbinary, 
#        width = width, height = height, dpi = dpi)
# 
# ggsave("hypothesis1_updated_RR_no_controls.png", 
#        plot = change_in_coefficients_plot_nonbinary_no_controls, 
#        width = width, height = height, dpi = dpi)

# Index outcomes, country image------------------------------------------------------
# Define your models
models_france2 <- list(
  mod1 = lm(image_chn_index ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_france),
  mod2 = lm(fp_cooperation_chn_index ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_france),
  mod3 = lm(image_usa_index ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_france),
  mod4 = lm(fp_cooperation_usa_index ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_france)
)

models_france2_no_controls <- list(
  mod1 = lm(image_chn_index ~ TREATMENT_relevel2, data = dt_chn_france),
  mod2 = lm(fp_cooperation_chn_index ~ TREATMENT_relevel2, data = dt_chn_france),
  mod3 = lm(image_usa_index ~ TREATMENT_relevel2, data = dt_chn_france),
  mod4 = lm(fp_cooperation_usa_index ~ TREATMENT_relevel2, data = dt_chn_france)
)

# Define your models
models_brazil2 <- list(
  mod1 = lm(image_chn_index ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_brazil),
  mod2 = lm(fp_cooperation_chn_index ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_brazil),
  mod3 = lm(image_usa_index ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_brazil),
  mod4 = lm(fp_cooperation_usa_index ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_brazil)
)

models_brazil2_no_controls <- list(
  mod1 = lm(image_chn_index ~ TREATMENT_relevel2, data = dt_chn_brazil),
  mod2 = lm(fp_cooperation_chn_index ~ TREATMENT_relevel2, data = dt_chn_brazil),
  mod3 = lm(image_usa_index ~ TREATMENT_relevel2, data = dt_chn_brazil),
  mod4 = lm(fp_cooperation_usa_index ~ TREATMENT_relevel2, data = dt_chn_brazil)
)

# Generate model plot data for France
nb_france2 <- modelplot(models_france2, coef_map = cm, draw = FALSE)
nb_france2_no_controls <- modelplot(models_france2_no_controls, coef_map = cm_no_controls, draw = FALSE)

# Generate model plot data for Brazil
nb_brazil2 <- modelplot(models_brazil2, coef_map = cm, draw = FALSE)
nb_brazil2_no_controls <- modelplot(models_brazil2_no_controls, coef_map = cm_no_controls, draw = FALSE)

# Combine the French and Brazilian data, adding a 'Country' column
nb_combined2 <- bind_rows(
  mutate(nb_france2, Country = "France"),
  mutate(nb_brazil2, Country = "Brazil")
)

nb_combined2_no_controls <- bind_rows(
  mutate(nb_france2_no_controls, Country = "France"),
  mutate(nb_brazil2_no_controls, Country = "Brazil")
)

# Assign Colors
new_colors <- c("Confidence" = "#8da0cb", "Trust" = "#fc8d62", 
                "Cooperation" = "#66c2a5", "Diplomacy" = "#e78ac3", 
                "Aid" = "#a6d854", "Business" = "#ffd92f", 
                "Leader" = "#e5c494")

# Renaming the models (adjust as necessary)
nb_combined2 <- nb_combined2 %>% 
  dplyr::mutate(model = dplyr::recode(model, 
                        "mod1" = "China Image", 
                        "mod2" = "China FP Cooperation",
                        "mod3" = "US Image", 
                        "mod4" = "US FP Cooperation"
  ))

nb_combined2_no_controls <- nb_combined2_no_controls %>% 
  dplyr::mutate(model = dplyr::recode(model, 
                                      "mod1" = "China Image", 
                                      "mod2" = "China FP Cooperation",
                                      "mod3" = "US Image", 
                                      "mod4" = "US FP Cooperation"
  ))


nb_combined2 <- nb_combined2 %>% 
  filter(term == "China UN Leadership" | term == "US UN Leadership")

nb_combined2_no_controls <- nb_combined2_no_controls %>% 
  filter(term == "China UN Leadership" | term == "US UN Leadership")

# Reorder the Country factor to ensure France is first
nb_combined2$Country <- factor(nb_combined2$Country, levels = c("France", "Brazil"))
nb_combined2_no_controls$Country <- factor(nb_combined2_no_controls$Country, levels = c("France", "Brazil"))

new_colors2 <- c("China Image" = "#8da0cb", "US Image" = "#fc8d62", 
                 "China FP Cooperation" = "#66c2a5", "US FP Cooperation" = "#e78ac3")

# Create the plot with facet wrap for each country
change_in_coefficients_plot_nonbinary2 <- ggplot(nb_combined2, 
       aes(y = term, x = estimate, 
       xmin = conf.low, xmax = conf.high, 
       color = model)) +
  geom_pointrange(position = position_dodge2(width = 0.5, reverse = TRUE), 
                  size = 0.25) +
  geom_vline(xintercept = 0) +
  scale_color_manual(values = new_colors2) +
  labs(x = "", y = "", 
       title = "",
       caption = "Bars represent 95% confidence intervals for coefficient estimates",
       color = "Model") +
  facet_wrap(~Country, ncol = 1) +
  theme_bw() +
  theme(#plot.title = element_text(hjust = 0.5, family = "Times New Roman"),
        plot.caption = element_text(hjust = 0, family = "Times New Roman"),
        legend.position = "bottom",
        legend.title = element_text(family = "Times New Roman"),
        legend.text = element_text(family = "Times New Roman"),
        axis.title = element_text(family = "Times New Roman"),
        axis.text = element_text(family = "Times New Roman", size = 10))

# Figure 2
change_in_coefficients_plot_nonbinary2

change_in_coefficients_plot_nonbinary2_no_controls <- ggplot(nb_combined2_no_controls, 
                                                 aes(y = term, x = estimate, 
                                                     xmin = conf.low, xmax = conf.high, 
                                                     color = model)) +
  geom_pointrange(position = position_dodge2(width = 0.5, reverse = TRUE), 
                  size = 0.25) +
  geom_vline(xintercept = 0) +
  scale_color_manual(values = new_colors2) +
  labs(x = "", y = "", 
       title = "",
       caption = "Bars represent 95% confidence intervals for coefficient estimates",
       color = "Model") +
  facet_wrap(~Country, ncol = 1) +
  theme_bw() +
  theme(#plot.title = element_text(hjust = 0.5, family = "Times New Roman"),
    plot.caption = element_text(hjust = 0, family = "Times New Roman"),
    legend.position = "bottom",
    legend.title = element_text(family = "Times New Roman"),
    legend.text = element_text(family = "Times New Roman"),
    axis.title = element_text(family = "Times New Roman"),
    axis.text = element_text(family = "Times New Roman", size = 10))

# Figure A-5
change_in_coefficients_plot_nonbinary2_no_controls

# Save the plot using ggsave for high resolution
# width <- 8 # in inches
# height <- 6 # in inches
# dpi <- 300 # adjust for desired resolution
# ggsave("hypothesis1_index_outcomes.png", 
#        plot = change_in_coefficients_plot_nonbinary2, 
#        width = width, height = height, dpi = dpi)
# 
# ggsave("hypothesis1_index_outcomes_no_controls.png", 
#        plot = change_in_coefficients_plot_nonbinary2_no_controls, 
#        width = width, height = height, dpi = dpi)

# Manipulation check -----------------------------------------------------
#manip check
dt_chn_france_manip <- dt_chn_france %>%
  mutate(.,
         #manip: 
         manip_china = case_when(
           MANIP_CHECK=="China"~1,
           MANIP_CHECK=="Germany"~0,
           MANIP_CHECK=="A Different Country"~0,
           MANIP_CHECK=="Not Mentioned"~0,
           MANIP_CHECK=="Switzerland"~0,
           MANIP_CHECK=="The United States"~0),
         
         manip_usa = case_when(
           MANIP_CHECK=="China"~0,
           MANIP_CHECK=="Germany"~0,
           MANIP_CHECK=="A Different Country"~0,
           MANIP_CHECK=="Not Mentioned"~0,
           MANIP_CHECK=="Switzerland"~0,
           MANIP_CHECK=="The United States"~1),
         
         #treatment
         treat_china = case_when(
           TREATMENT_relevel2=="1"~1,
           TREATMENT_relevel2=="2"~0,
           TREATMENT_relevel2=="Control"~0),
         
         treat_usa = case_when(
           TREATMENT_relevel2=="1"~0,
           TREATMENT_relevel2=="2"~1,
           TREATMENT_relevel2=="Control"~0),
         
         treat_control = case_when(
           TREATMENT_relevel2=="1"~0,
           TREATMENT_relevel2=="2"~1,
           TREATMENT_relevel2=="Control"~1)
         )

man<- lm(manip_china~treat_china,dt_chn_france_manip)
summary(man) #resp in the china condition more likely to (correctly) say head was china

man2<- lm(manip_usa~treat_usa,dt_chn_france_manip)
summary(man2) #resp in the usa condition more likely to (correctly) say head was usa

dt_chn_brazil_manip <- dt_chn_brazil %>%
  mutate(.,
         #manip: 
         manip_china = case_when(
           MANIP_CHECK=="China"~1,
           MANIP_CHECK=="Germany"~0,
           MANIP_CHECK=="A Different Country"~0,
           MANIP_CHECK=="Not Mentioned"~0,
           MANIP_CHECK=="Switzerland"~0,
           MANIP_CHECK=="The United States"~0),
         
         manip_usa = case_when(
           MANIP_CHECK=="China"~0,
           MANIP_CHECK=="Germany"~0,
           MANIP_CHECK=="A Different Country"~0,
           MANIP_CHECK=="Not Mentioned"~0,
           MANIP_CHECK=="Switzerland"~0,
           MANIP_CHECK=="The United States"~1),
         
         #treatment
         treat_china = case_when(
           TREATMENT_relevel2=="1"~1,
           TREATMENT_relevel2=="2"~0,
           TREATMENT_relevel2=="Control"~0),
         
         treat_usa = case_when(
           TREATMENT_relevel2=="1"~0,
           TREATMENT_relevel2=="2"~1,
           TREATMENT_relevel2=="Control"~0),
         
         treat_control = case_when(
           TREATMENT_relevel2=="1"~0,
           TREATMENT_relevel2=="2"~0,
           TREATMENT_relevel2=="Control"~1)
  )

man<- lm(manip_china~treat_china,dt_chn_brazil_manip)
summary(man) #resp in the china condition more likely to (correctly) say head was china

man2<- lm(manip_usa~treat_usa,dt_chn_brazil_manip)
summary(man2) #resp in the usa condition more likely to (correctly) say head was usa

# Test: compare to 'no country' baseline for state image----------------------------------------------
# Create new treatment indicator using manip check
test_data_france <- dt_chn_france %>% 
  filter(MANIP_CHECK %in% "Not Mentioned" |
           MANIP_CHECK %in% "China" |
           MANIP_CHECK %in% "The United States")

test_data_france$TREATMENT_manip <- relevel(as.factor(test_data_france$MANIP_CHECK), 
                                            ref = "Not Mentioned")

test_data_brazil <- dt_chn_brazil %>% 
  filter(MANIP_CHECK %in% "Not Mentioned" |
           MANIP_CHECK %in% "China" |
           MANIP_CHECK %in% "The United States")

test_data_brazil$TREATMENT_manip <- relevel(as.factor(test_data_brazil$MANIP_CHECK), 
                                            ref = "Not Mentioned")

# Model Summary UN Legitimacy-----------------------------------------------------
# Define your models
models <- list(
  mod1 = lm(un_index ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt),
  mod2 = lm(un_index ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_france),
  mod3 = lm(un_index ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_brazil)
)

models_no_controls <- list(
  mod1 = lm(un_index ~ TREATMENT_relevel2, data = dt),
  mod2 = lm(un_index ~ TREATMENT_relevel2, data = dt_chn_france),
  mod3 = lm(un_index ~ TREATMENT_relevel2, data = dt_chn_brazil)
)

# Generate model plot data
nb <- modelplot(models, coef_map = cm, draw = FALSE)
nb_no_controls <- modelplot(models_no_controls, coef_map = cm_no_controls, draw = FALSE)

# Renaming the models (adjust as necessary)
nb <- nb %>% 
  dplyr::mutate(model = dplyr::recode(model, 
                        "mod1" = "Overall", 
                        "mod2" = "France",
                        "mod3" = "Brazil", 
  ))

nb_no_controls <- nb_no_controls %>% 
  dplyr::mutate(model = dplyr::recode(model, 
                                      "mod1" = "Overall", 
                                      "mod2" = "France",
                                      "mod3" = "Brazil", 
  ))


# Assign Colors
new_colors <- c("Overall" = "#8da0cb", "France" = "#fc8d62", 
                "Brazil" = "#66c2a5")

nb <- nb %>% filter(term == "China UN Leadership" | term == "US UN Leadership")
nb_no_controls <- nb_no_controls %>% filter(term == "China UN Leadership" | term == "US UN Leadership")

# Create the plot with facet wrap for each country
change_in_coefficients_plot_nonbinary <- ggplot(nb, aes(y = term, x = estimate, 
                                                                 xmin = conf.low, xmax = conf.high, 
                                                                 color = model)) +
  geom_pointrange(position = position_dodge2(width = 0.5, reverse = TRUE), size = 0.25) +
  geom_vline(xintercept = 0) +
  scale_color_manual(values = new_colors) +
  labs(x = "IO Legitimacy", y = "", 
       title = "Effect of China's UN Leadership on IO Legitimacy",
       caption = "Bars represent 95% confidence intervals for coefficient estimates",
       color = "Model") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, family = "Times New Roman"),
        plot.caption = element_text(hjust = 0, family = "Times New Roman"),
        legend.position = "bottom",
        legend.title = element_text(family = "Times New Roman"),
        legend.text = element_text(family = "Times New Roman"),
        axis.title = element_text(family = "Times New Roman"),
        axis.text = element_text(family = "Times New Roman", size = 10))

# Figure 3
change_in_coefficients_plot_nonbinary

change_in_coefficients_plot_nonbinary_no_controls <- ggplot(nb_no_controls, aes(y = term, x = estimate, 
                                                        xmin = conf.low, xmax = conf.high, 
                                                        color = model)) +
  geom_pointrange(position = position_dodge2(width = 0.5, reverse = TRUE), size = 0.25) +
  geom_vline(xintercept = 0) +
  scale_color_manual(values = new_colors) +
  labs(x = "IO Legitimacy", y = "", 
       title = "Effect of China's UN Leadership on IO Legitimacy",
       caption = "Bars represent 95% confidence intervals for coefficient estimates",
       color = "Model") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, family = "Times New Roman"),
        plot.caption = element_text(hjust = 0, family = "Times New Roman"),
        legend.position = "bottom",
        legend.title = element_text(family = "Times New Roman"),
        legend.text = element_text(family = "Times New Roman"),
        axis.title = element_text(family = "Times New Roman"),
        axis.text = element_text(family = "Times New Roman", size = 10))

# Figure A-7
change_in_coefficients_plot_nonbinary

# Save the plot using ggsave
# ggsave("hypothesis2_updated.png", plot = change_in_coefficients_plot_nonbinary, 
#        width = width, height = height, dpi = dpi)
# 
# ggsave("hypothesis2_updated_no_controls.png", 
#        plot = change_in_coefficients_plot_nonbinary_no_controls, width = width, height = height, dpi = dpi)

# Individual UN outcomes --------------------------------------------
# Define your models
models_un <- list(
  mod1 = lm(leg_un ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_france),
  mod2 = lm(trust_un ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_france),
  mod3 = lm(rep1_un ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_france),
  mod4 = lm(rep2_un ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_france)
)

models_un_no_controls <- list(
  mod1 = lm(leg_un ~ TREATMENT_relevel2, data = dt_chn_france),
  mod2 = lm(trust_un ~ TREATMENT_relevel2, data = dt_chn_france),
  mod3 = lm(rep1_un ~ TREATMENT_relevel2, data = dt_chn_france),
  mod4 = lm(rep2_un ~ TREATMENT_relevel2, data = dt_chn_france)
)

# Define your models
models_un2 <- list(
  mod1 = lm(leg_un ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_brazil),
  mod2 = lm(trust_un ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_brazil),
  mod3 = lm(rep1_un ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_brazil),
  mod4 = lm(rep2_un ~ TREATMENT_relevel2 + education + 
              income + 
              Age + male +
              trustingov + polinterest +
              fp_index + ideol +
              chn_friendly + chn_threat +
              usa_friendly + usa_threat, data = dt_chn_brazil)
)

models_un2_no_controls <- list(
  mod1 = lm(leg_un ~ TREATMENT_relevel2, data = dt_chn_brazil),
  mod2 = lm(trust_un ~ TREATMENT_relevel2, data = dt_chn_brazil),
  mod3 = lm(rep1_un ~ TREATMENT_relevel2, data = dt_chn_brazil),
  mod4 = lm(rep2_un ~ TREATMENT_relevel2, data = dt_chn_brazil)
)

# Generate model plot data for France
nb_un2 <- modelplot(models_un, coef_map = cm, draw = FALSE)
nb_un2_no_controls <- modelplot(models_un_no_controls, coef_map = cm_no_controls, draw = FALSE)

# Generate model plot data for Brazil
nb_brazil2 <- modelplot(models_un2, coef_map = cm, draw = FALSE)
nb_brazil2_no_controls <- modelplot(models_un2_no_controls, coef_map = cm_no_controls, draw = FALSE)

# Combine the French and Brazilian data, adding a 'Country' column
nb_combined2 <- bind_rows(
  dplyr::mutate(nb_un2, Country = "France"),
  dplyr::mutate(nb_brazil2, Country = "Brazil")
)

nb_combined2_no_controls <- bind_rows(
  dplyr::mutate(nb_un2_no_controls, Country = "France"),
  dplyr::mutate(nb_brazil2_no_controls, Country = "Brazil")
)

nb_combined2 <- nb_combined2 %>% 
  dplyr::mutate(model = dplyr::recode(model, 
                        "mod1" = "Confidence", 
                        "mod2" = "Trust",
                        "mod3" = "Reputation",
                        "mod4" = "Approve"
  ))

nb_combined2_no_controls <- nb_combined2_no_controls %>% 
  dplyr::mutate(model = dplyr::recode(model, 
                                      "mod1" = "Confidence", 
                                      "mod2" = "Trust",
                                      "mod3" = "Reputation",
                                      "mod4" = "Approve"
  ))

# Assign Colors
new_colors <- c("Confidence" = "#8da0cb", "Trust" = "#fc8d62", 
                "Reputation" = "#66c2a5", "Approve" = "#e78ac3")

nb_combined2 <- nb_combined2 %>% 
  filter(term == "China UN Leadership" | term == "US UN Leadership")

nb_combined2_no_controls <- nb_combined2_no_controls %>% 
  filter(term == "China UN Leadership" | term == "US UN Leadership")
# Reorder the Country factor to ensure France is first
nb_combined2$Country <- factor(nb_combined2$Country, levels = c("France", "Brazil"))

nb_combined2_no_controls$Country <- factor(nb_combined2_no_controls$Country, levels = c("France", "Brazil"))

# Create the plot with facet wrap for each country
change_in_coefficients_plot_nonbinary2 <- ggplot(nb_combined2, 
                                                 aes(y = term, x = estimate, 
                                                     xmin = conf.low, xmax = conf.high, 
                                                     color = model)) +
  geom_pointrange(position = position_dodge2(width = 0.5, reverse = TRUE), 
                  size = 0.25) +
  geom_vline(xintercept = 0) +
  scale_color_manual(values = new_colors) +
  labs(x = "", y = "", 
       title = "",
       caption = "Bars represent 95% confidence intervals for coefficient estimates",
       color = "Model") +
  facet_wrap(~Country, ncol = 1) +
  theme_bw() +
  theme(#plot.title = element_text(hjust = 0.5, family = "Times New Roman"),
    plot.caption = element_text(hjust = 0, family = "Times New Roman"),
    legend.position = "bottom",
    legend.title = element_text(family = "Times New Roman"),
    legend.text = element_text(family = "Times New Roman"),
    axis.title = element_text(family = "Times New Roman"),
    axis.text = element_text(family = "Times New Roman", size = 10))

# Figure A-4
change_in_coefficients_plot_nonbinary2

change_in_coefficients_plot_nonbinary2_no_controls <- ggplot(nb_combined2_no_controls, 
                                                 aes(y = term, x = estimate, 
                                                     xmin = conf.low, xmax = conf.high, 
                                                     color = model)) +
  geom_pointrange(position = position_dodge2(width = 0.5, reverse = TRUE), 
                  size = 0.25) +
  geom_vline(xintercept = 0) +
  scale_color_manual(values = new_colors) +
  labs(x = "", y = "", 
       title = "",
       caption = "Bars represent 95% confidence intervals for coefficient estimates",
       color = "Model") +
  facet_wrap(~Country, ncol = 1) +
  theme_bw() +
  theme(#plot.title = element_text(hjust = 0.5, family = "Times New Roman"),
    plot.caption = element_text(hjust = 0, family = "Times New Roman"),
    legend.position = "bottom",
    legend.title = element_text(family = "Times New Roman"),
    legend.text = element_text(family = "Times New Roman"),
    axis.title = element_text(family = "Times New Roman"),
    axis.text = element_text(family = "Times New Roman", size = 10))

# Figure A-8
change_in_coefficients_plot_nonbinary2

# Save the plot using ggsave for high resolution
# ggsave("hypothesis2_indivdual_outcomesRR.png", 
#        plot = change_in_coefficients_plot_nonbinary2, 
#        width = width, height = height, dpi = dpi)
# 
# ggsave("hypothesis2_indivdual_outcomesRR_no_controls.png", 
#        plot = change_in_coefficients_plot_nonbinary2_no_controls, 
#        width = width, height = height, dpi = dpi)

# Analysis of in-text refs to Parizek data  --------------------------------------------
parizek <- read.csv("IGOs_Media_Visibility_LocalLanguages_diso.csv")

# attention to UN
parizek_trim <- parizek %>% 
  dplyr::select(country_name_ggplot_map, UNO) 

parizek_trim <- parizek_trim %>% 
  arrange(desc(UNO)) %>%
  mutate(rank = 1:nrow(parizek),
         percent_rank = percent_rank(UNO))

parizek_trim$UNO[parizek_trim$country_name_ggplot_map %in% "Brazil"]
parizek_trim$rank[parizek_trim$country_name_ggplot_map %in% "Brazil"]
parizek_trim$percent_rank[parizek_trim$country_name_ggplot_map %in% "Brazil"]

parizek_trim$UNO[parizek_trim$country_name_ggplot_map %in% "France"]
parizek_trim$rank[parizek_trim$country_name_ggplot_map %in% "France"]
parizek_trim$percent_rank[parizek_trim$country_name_ggplot_map %in% "France"]

summary(parizek_trim$UNO)
quantile(parizek_trim$UNO)

# compare to attention to all IOs
parizek_trim_all <- parizek %>% 
  dplyr::select(country_name_ggplot_map, igos_all) 

parizek_trim_all <- parizek_trim_all %>% 
  arrange(desc(igos_all)) %>%
  mutate(rank = 1:nrow(parizek),
         percent_rank = percent_rank(igos_all))

parizek_trim_all$igos_all[parizek_trim_all$country_name_ggplot_map %in% "Brazil"]
parizek_trim_all$rank[parizek_trim_all$country_name_ggplot_map %in% "Brazil"]
parizek_trim_all$percent_rank[parizek_trim_all$country_name_ggplot_map %in% "Brazil"]

parizek_trim_all$igos_all[parizek_trim_all$country_name_ggplot_map %in% "France"]
parizek_trim_all$rank[parizek_trim_all$country_name_ggplot_map %in% "France"]
parizek_trim_all$percent_rank[parizek_trim_all$country_name_ggplot_map %in% "France"]

summary(parizek_trim_all$igos_all)
quantile(parizek_trim_all$igos_all)


# Figure A-1
# Import pre-cleaned data #from Mattingly et al https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/CQ4FZR
dat_cn <- import("cgtn.csv", encoding = "UTF-8") 

# Convert the data to character format for string operations
dat_cn_text <- apply(dat_cn, 2, as.character)

# Collapse all text columns into a single character vector
dat_cn_combined <- paste(dat_cn_text, collapse = " ")

# Convert the data to character format for string operations
dat_cn_text <- apply(dat_cn, 2, as.character)

# Collapse all text columns into a single character vector
dat_cn_combined <- paste(dat_cn_text, collapse = " ")

# Define a vector of keywords related to the United Nations and other soft power efforts
un_keywords <- c("olympics", "tianwen", "expo", "UN", "WHO", "UNDP")

# Count occurrences for each keyword
keyword_counts <- sapply(un_keywords, function(keyword) 
  stringr::str_count(dat_cn_combined, paste0("\\b", keyword, "\\b")))

# Convert to a data frame
keyword_df <- data.frame(Keyword = un_keywords, Count = keyword_counts)

# Rename "tianwen" to "space"
keyword_df$Keyword <- ifelse(keyword_df$Keyword == "tianwen", "Space", keyword_df$Keyword)
keyword_df$Keyword <- ifelse(keyword_df$Keyword == "expo", "Expo", keyword_df$Keyword)
keyword_df$Keyword <- ifelse(keyword_df$Keyword == "olympics", "Olympics", keyword_df$Keyword)

# Define UN-related keywords
un_related <- c("UN", "WHO", "UNDP")

# Assign colors: UN organizations in UN Blue, others in Gray
keyword_df$Color <- ifelse(keyword_df$Keyword %in% un_related, "#009EDB", "#A9A9A9")

# Create the bar plot with custom font
output_file <- ggplot(keyword_df, aes(x = reorder(Keyword, -Count), y = Count, fill = Keyword)) +
  geom_bar(stat = "identity", show.legend = FALSE) +
  scale_fill_manual(values = setNames(keyword_df$Color, keyword_df$Keyword)) +  # Apply color mapping
  labs(title = "China's CGTV Global Coverage",
       x = "Keyword",
       y = "Count") +
  theme_minimal() +
  theme(
    text = element_text(family = "Times New Roman"), 
    axis.text.x = element_text(angle = 45, hjust = 1, family = "Times New Roman"),  
    axis.text.y = element_text(family = "Times New Roman"),  
    title = element_text(family = "Times New Roman")  
  )

output_file
# # Set dimensions and resolution (adjust as needed)
# width <- 8 # in inches
# height <- 6 # in inches
# dpi <- 300 # adjust for desired resolution
# 
# 
# # Save the plot using ggsave for high resolution
# ggsave("descriptive.png", plot = output_file, width = width, height = height, dpi = dpi)
# 


# Figure A-2
# create a data frame with the leader tenure information
leaders <- data.frame(
  Executive_Head = c("Carlos Magariños", "Kandeh Yumkella", "Yong Li", "Gerd Müller"),
  Nationality = c("Argentina", "Sierra Leone", "China", "Germany"),
  Term_Start = c(1997, 2005, 2013, 2021),
  Term_End = c(2005, 2012, 2021, 2023)
)

# plot the leaders' tenure using ggplot
UNIDO <- ggplot(leaders, aes(x = Term_Start, xend = Term_End, y = Executive_Head, yend = Executive_Head)) +
  geom_segment(aes(color = Nationality), size = 3) +
  geom_point(aes(color = Nationality), size = 5, shape = 21, fill = "white") +
  scale_color_manual(values = c("Argentina" = "#6C3082", "Sierra Leone" = "#008080", 
                                "China" = "red", "Germany" = "#009E73")) +
  labs(x = "Year", y = "", title = "United Nations Industrial Development Organization") +
  theme_minimal() + theme(text = element_text(size=16)) 


#FAO
# create a data frame with the leader tenure information
fao_leaders <- data.frame(
  Executive_Head = c("Jacques Diouf", "José Graziano da Silva", "Dongyu Qu"),
  Nationality = c("Senegal", "Brazil", "China"),
  Term_Start = c(1994, 2012, 2019),
  Term_End = c(2011, 2019, 2023)
)

FAO <- ggplot(fao_leaders, aes(x = Term_Start, xend = Term_End, y = Executive_Head, yend = Executive_Head)) +
  geom_segment(aes(color = Nationality), size = 3) +
  geom_point(aes(color = Nationality), size = 5, shape = 21, fill = "white") +
  scale_color_manual(values = c("Senegal" = "#6C3082", "Brazil" = "#008080", 
                                "China" = "red")) +
  labs(x = "Year", y = "", title = "Food and Agriculture Organization") +
  theme_minimal() + theme(text = element_text(size=16)) 

itu_leaders <- data.frame(
  Executive_Head = c("Yoshio Utsumi", "Dr Hamadoun Touré", "Houlin Zhao", "Doreen Bogden Martin"),
  Nationality = c("Japan", "Mali", "China", "United States"),
  Term_Start = c(1998, 2007, 2015, 2023),
  Term_End = c(2002, 2014, 2022, 2023)
)

ITU <- ggplot(itu_leaders, aes(x = Term_Start, xend = Term_End, y = Executive_Head, yend = Executive_Head)) +
  geom_segment(aes(color = Nationality), size = 3) +
  geom_point(aes(color = Nationality), size = 5, shape = 21, fill = "white") +
  scale_color_manual(values = c("Japan" = "#6C3082", "Mali" = "#008080", 
                                "China" = "red", "United States" = "blue")) +
  labs(x = "Year", y = "", title = "International Telecommunications Union") +
  theme_minimal() + theme(text = element_text(size=16)) 

#ICAO
icao_leaders <- data.frame(
  Executive_Head = c("Taïeb Chérif", "Raymond Benjamin", "Dr. Fang Liu", "Juan Carlos Salazar"),
  Nationality = c("Algeria", "France", "China", "Colombia"),
  Term_Start = c(2003, 2009, 2015, 2021),
  Term_End = c(2009, 2015, 2021, 2023)
)

ICAO <- ggplot(icao_leaders, aes(x = Term_Start, xend = Term_End, y = Executive_Head, yend = Executive_Head)) +
  geom_segment(aes(color = Nationality), size = 3) +
  geom_point(aes(color = Nationality), size = 5, shape = 21, fill = "white") +
  scale_color_manual(values = c("Algeria" = "#6C3082", "France" = "#008080", 
                                "China" = "red", "Colombia" = "#009E73")) +
  labs(x = "Year", y = "", title = "International Civil Aviation Organization") +
  theme_minimal() + theme(text = element_text(size=16)) 

ggarrange(UNIDO, FAO, ITU, ICAO,
          ncol = 2, nrow = 2)
